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We present an analytical scheme, easily implemented nu- 
merically, to generate synthetic Gaussian turbulent flows by 
using a linear Langevin equation, where the noise term acts 
as a stochastic stirring force. The characteristic parameters 
of the velocity field are well introduced, in particular the kine- 
matic viscosity and the spectrum of energy. As an application, 
the diffusion of a passive scalar is studied for two different en- 
ergy spectra. Numerical results are compared favorably with 
analytical calculations. 



PACS numbers: 05.40+j, 47.27.Eq, 47.27.Qb 



I. INTRODUCTION 

The statistical approach to turbulence has a long his- 
tory on its own (!],§]. Actually, it is the firm recogni- 
tion that both fundamental and applied aspects of turbu- 
lence can be consistently addressed from the knowledge 
of the statistical moments and correlations of the velocity 
flow that has made this approach particularly appealing. 
Since in any case such an statistical point of view has 
adopted multiple perspectives during these last decades, 
let us first put our formulation into an appropriate con- 
text. 

Our primary aim here is to describe an analytical 
and numerical methodology to generate a statistically 
stationary, homogeneous, and isotropic two-dimensional 
turbulent flow with zero mean velocity and well-defined 
energy spectrum. Needless to say, such a regime of steady 
turbulence can only be maintained by means of an exter- 
nal input of energy to compensate the dissipative na- 
ture of the viscous forces. This external input could 
be conceptually associated with stirring forces which are 
stochastic in nature so as to produce a random velocity 
field which is to represent the turbulent flow |J. It is 
quite obvious that the statistical properties of the steady 
stirring forces will ultimately determine those of the tur- 
bulent flow. 

Ideally, the statistical properties of any turbulent flow 
should come as the output of a first principles, Navier- 
Stokes based, formulation of the problem. However, we 
will adopt here a somewhat reversed perspective aimed 
at developing a methodology to construct what would 



be a sort of synthetic turbulence [|4|-|3|. Rather than to 
retain the nonlinear coupling which makes possible the 
redistribution of energy from the largest length scales 
down into the smaller ones, we assume that the energy 
is incorporated into the system in an individual wave 
number basis. In other words, our model would represent 
a collection of uncoupled stirrers each one acting on a 
single-length scale and introducing its own wave number 
dependent energy contribution, but chosen to reproduce, 
and this is the main point in our approach, the desired 
spectral distribution in steady state. 

This practical approach to produce a turbulent flow 
is justified by the fact that what we want is to study a 
physical process within the turbulent medium described 
only by its statistical properties. Our strategy is accom- 
plished by using a generalized Ornstein Uhlenbeck-like 
Langevin equation for the stream function r?(r,t) |||3j. 
Being more specific, our proposal is based on the use of 
the following Langevin equation 



dt 



= ^V 2 77(r,t) + Q[A 2 V 2 ]V.C(r,i), 



(1) 



where v is the kinematic viscosity and C(r, t) is a Gaus- 
sian white noise with zero mean value and correlation 

(C'(ri,ii)C J '(r 2 ,t 2 )) =2 e o^(ti-i2^(ri-r 2 )^'. (2) 

In this last expression, eo, and A are parameters which 
control, respectively, the intensity and correlation length 
of the random flow. Q[A 2 V 2 ] will play the most rele- 
vant role in our scheme, representing the random stirring 
forces. 

The random flow generated in this way has Gaussian 
properties due to the linear nature of such an equation. 
Apart from that, other limitations concerning some in- 
variances of the turbulent flow appear also as a conse- 
quence of such a linearity ||. However, the big advan- 
tage of our model is that it facilitates the control over 
the characteristic parameters of the turbulence, i.e. its 
integral time and length scales and its spectrum, just 
by appropriately prescribing the input parameters of the 
noise entering into the Langevin equation. A final remark 
concerning Gaussianity is worth mentioning at this point. 
Certainly, we recognize that such an statistical property 
is under scrutiny j9|. However particular experimental 
scenarios supporting it flO[| , together with the lacking of 
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conclusive results on intermittency 0, leave this ques- 
tion rather open and allows us to use such a Gaussian 
property, at least as a working simplifying hypothesis, 
mainly when the focus is in the study of physical process 
inside this medium. 

This paper is organized as follows. In the next section 
we summarize the Langevin approach to generate the 



random flow (synthetic turbulence). In Sec. Ill the tech- 
nical details of the numerical simulations are presented. 
As a sort of application we consider the classical problem 



of the diffusion of a passive scalar in Sec. IV. We devote 
Sec. |V| to draw some conclusions and perspectives. 

II. LANGEVIN APPROACH 

As our approach starts with the generation of a scalar 
stream function jllj , we are going to review first the re- 
lationships between the stream function properties and 
those of the velocity field. Let ry(r, t) be the stream func- 
tion from which we define the two-dimensional incom- 
pressible turbulent field v(r, t) (V.v = 0), 



v(r,i) 



dr](r,t) dr](r,t) 



dy 



dx 



(3) 



Its mean value is to be taken zero and the two-point 
velocity autocorrelation is defined as usual through 

(« i (ri,tiy(r 2 ,t 2 ))=i2«'(ri-r 2 ,t 1 -t2) ) (4) 

where homogeneity in space and time is explicitly de- 
noted |^,^|. Our scheme is entirely implemented in two 
dimensions (2D). However, it could be generalized to 3D 
just by taking a vector stream function T] and obtaining 
the velocity field as v = V x rj. In this case each com- 
ponent of the stream function must satisfy a Langevin 
equation similar to (Q) with independent noises. 

Making use of the spatial isotropy we define the radial 
correlation function as 



R(r,s) = -[R **{r,s)+R«y(r,s)}, 



(5) 



where r = |ri — r 2 1 and s = \tx — 1 2 |- The autocorrelation 
of T?(r, t) 



C(r, s) = (77(ri,ii)77(r 2 ,i 2 )) 



(6) 



is correspondingly assumed to have the properties of ho- 
mogeneity, isotropy, and stationarity. From the definition 
(|J) we can express the velocity correlation (||) in terms 

of M) 



R(r,s) = 



4tt 



dkk 3 J (kr)C(k,s), 



(7) 



where Jo(kr) is the Bessel function of zeroth order and 
C(k, s) is the Fourier transform of C(r, s). 



The physical parameters of steady turbulence, i.e., 
its intensity it 2 , and characteristic (integral) time and 
length scales follow from their standard definitions pl[| , 



u 2 = R(0,0), 



1 f°° 
to = — j / ds i?(0, s), 
u o Jo 

dr R(r,0). 



(8) 



lo = 



x -'O 



With these definitions in mind let us move to the dis- 
cussion of the analytical scheme. The Fourier transform 
of (TO) reads 



dr](k,t) 
dt 



-vk 2 r)(k,t)-iQ{-\ 2 k 2 )k j C j (k,t). (9) 



Now we introduce the structure function S(k,t) related 
with the equal-time correlation function (^), but in 
Fourier space 

faCki, t)»7(k a , *)) - (27r) 2 5(ki + k 2 )5(fc, t). (10) 

By using standard stochastic calculus [Q, S(k,t) is 
shown to verify 



dS(k,t) 
dt 



-2vkfS(k, t) + 2e vk z Q z {-Yk z ). (11) 



On the other hand, C(k,s) obeys, in the steady state, 
the equation 



dC(k,s) 
ds 



-vk 2 C(k,s), 



(12) 



with the initial condition 

C(k,0) = S*(fc,t-> oo) = S st (k). (13) 
From (|ll|) and Q we get 

C(M = 6 Q 2 (-A 2 fc 2 )e- fe2s , (14) 

The energy spectrum E(k,t) is defined in terms of 
S(k,t) as 



E(k,t) = —k 3 S(k,t). 
in 



(15) 



Using (|T^) we can also obtain the equation of evolution 
ofE(kJ) [HJl| 



dE(k,t) 
dt 



-2uk 2 E(k,t) + 2k 2 W(k), (16) 



where 



W ( k ) = e -^k 3 Q 2 (~\ 2 k 2 ). 



(17) 



This quantity can be regarded as the input of energy due 
to the stirring forces. The stationary state is achieved 
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when the input term is exactly balanced by the dissipa- 
tion term, 



t+At 



dt'C» v (t')exp[(t + At-t')v CfiV ], 

(22) 



E(k, t^oo) = E(k) = -W(k) = ^k 3 Q 2 (-X 2 k 2 ) = ^k 3 C(k, 0). 

where Q^ v in the last two expressions denote the dis- 



(18) 



crete Fourier transform of the operator Q[X V ]. In the 



The stirring intensity u 2 can also be related with the Fourier space the derivative operators have been trans- 
spectrum l&ted into 



dk E(k), 



(19) 



and the characteristic time to (|9|), and the characteristic 
length (||) can be related as well with the viscosity v 
and A. Explicit relations will be obtained later on for 
particular spectra. 

As it is quite obvious, our scheme does not incorporate 
the non-linear term of the Navier-Stokes equation. No- 
tice, nevertheless, that this scheme is versatile enough to 
reproduce a large variety of energy spectra, just by ap- 
propriately prescribing the differential operator Q[A 2 V 2 ] 
in (Q) . In our approach the noisy term represents the stir- 
ring mechanism which feeds energy continuously into the 
system according to the chosen spectrum. At the same 
time this energy is dissipated at time and spatial scales 
controlled by v. As there are only linear terms in the 
Langevin equation, no kind of energy cascade is possible. 
The noise forces introduce the whole hierarchy of turbu- 
lent structures evolving according to their own time and 
length scales, but without any interaction between them. 



III. SIMULATION ALGORITHM 

We have chosen for the discretization of the real space 
a standard two-dimensional square lattice N x N with 
elementary unit spacing A in both directions. In most of 
our simulations N ~ 128 and A = 0.5, unless other val- 
ues are specified. The discrete Fourier space is discretized 
accordingly k = (k x ,k y ) — 2n / N A(pt,v) Q (Greek in- 
dices are used in Fourier space on all that follows). We 
can take advantage of the fact that (^|), when written in 
the discrete Fourier space, transforms into a set of linear 
and not coupled ordinary differential equations. In these 
circumstances the exact integration between t and t + At, 
(At = 0.1t Q ) gives 

^(f + At) = exp(vc pv At) rip^t) + 0p V (t) + Jnv(t), 

(20) 

in terms of new random variables j3^ v it) and 7^(4) are 
defined according to 



t+At 



dt' CL(0 exp [(t + At- t>c p 



(21) 



2 

A 2 " 



/2ttu\ /2ttv\ 
cos — - + cos I - 



, (23) 



exp 



N J 



L , exp 1 



N 



(24) 



The correlation of the random variables can be ex- 
pressed as 

</£„(*)/W*)> = Q^Q P ad%d x p ; < QUt'KrW) > 

pt+At pt+At 

/ dt' / dt" exp [(t + At- t')vc^ v + vc pa (t + At- t") 
Jo Jo 

and similar equations for the 7^ (t) . Using (||) , the sym- 
metry properties of the Q operator Q^ v — Q-^ v = Q*^ v , 
and further integrating and using the expressions defined 
above for the discrete operators in Fourier space we fi- 
nally have 



, (2. 



(0%(t)P P *(t)) = 2e Q 5 w 5 va N 2 Ql v c-l [exp (2^A<) - 1] 

(26) 

We can now construct an explicit expression for f3fj, v (t) 
adapted to the result just obtained 



1 — cos 



f3^ v (t) = 2e N 2 Q 2 c-^ [exp (2^A<) - 1] 



(27) 



where a^ v are Gaussian random numbers, which satisfy 
(a*^ v (t)a pa (t)) = 8fj, p 8 va . We would proceed analogously 
for 7 jU «(t). The correlation of the stream function in the 
steady state is 



1/2 



(v; v V P a) st = Ql v e (NA) 2 S^S. v 



(28) 



When dealing with the diffusion of the passive scalars in 
Sec. |y| we will always start our simulations in a steady 
configuration of the random flow. According to the re- 
sult (^8j) this is easily accomplished by taking as initial 
condition 



V(°) = Q^l /2 (NA)a^ v (Q). 



(29) 
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At each time step and after having generated the 
stream function in Fourier space we proceed to antitrans- 
form and using it in relation with the appropriate dis- 
cretized forms for the velocity field. We skip some details 
to refer directly to the the discrete version of the energy 
spectrum which finally reads B 



E„ 



27VA 3 



(M 2 



,2\i, 



) 2 

■ fJ,V 



sin 



/ 2nn\ „ ( 2ttv 

— - + sm 2 

V N J \N 

(30) 



IV. DIFFUSION OF PASSIVE SCALARS 

Before we start the study of the scalar diffusion, we 
will discuss the selection of spectra. There are several 
possible choices and between all of them we have selected 
two spectra which are well behaved in all the range of 
wave numbers. 



(37) 



Our second choice is the Karman-Obukhov's (KO) 
spectrum which was introduced p6| to study Kolgomorov 
turbulence with a long " — 5/3" tail in the spectrum for 
large k. To this end we select the family of Karman- 
Obukhov spectra, whose general behavior is 



E(k) oc ¥• 



k 2 
K o. 



-(5+3n)/6 



For convenience we have taken n 
Q operator is then 



(38) 

3. The choice of the 



Q[A 2 V 2 



1 



A 2 V 2 ) 



-7/6 



(39) 



In this case we obtain the following results for the three 
basic parameters 



9eo 
32ttA 4 : 



(40) 



A. Kraichnan's and Karman-Obukhov spectra 

The spectrum introduced by Kraichnan (K) [[l5| de- 
scribes turbulent velocity fields with a widely distributed 
band of excitations and a peak centered at well-defined 
wave number ~ ko 



E(k) oc A; 3 exp 



k 2 



In this case, the operator Q has to be chosen according 
to 



Q[A 2 V 2 ] = exp 



A 2 V' 



(32) 



where A = k Q 1 . Using (|l4|) and substituting in (0), the 
velocity correlation function in steady state reads 



R(r,s) 



8tt(A 2 + us) 2 



1 - 



4(A 2 + us) 



exp 



4(A 2 + us) 



From it, we can obtain Ug, and the integral time and 
length scales according to 



8vrA 4 ' 



_ A 2 



(34) 
(35) 



A 2 

to — 7j— , 

3u 



(41) 



Zo = A 



r(i/2)r(5/6) 
2r(i/3) ' 



(42) 



(32) The corresponding Reynolds number is expressed as 



"-if 1 ' 



1 3r(5/6) 
4A^2r(l/3)' 



(43) 



A closed expression for the correlation function cannot 
be obtained in this case. 

These two spectra have been simulated according to 
the recipes of the previous section. In Fig. 1 we have 
plotted the evolution of an initially flat spectrum towards 
its final steady state. By looking at the patterns corre- 
sppjgMng to different time intervals we can judge the role 
of the kinematic viscosity. In particular the modes with 
smaller k relax slower than modes with larger k. Contin- 
uous and discrete analytical results are also compared in 
the steady state to see the differences introduced by the 
discretization procedure specially for large k. The clear 
differences between both spectra are going to influence 
the diffusion of passive scalar as it will be discussed later 



aVt? 



(36) 



The dimensionless Reynolds number, defined according 
to Re = IqUq/u, is thus expressed in terms of the noise 
parameters as 



B. Scalar diffusion in a continuous scheme 

Scalar diffusion in turbulent flows is a classical prob- 
lem |L7j which still deserves a lot of attention Our 
starting point is the Eulerian equation of motion for a 
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scalar distribution ip(r,t) which is assumed to be pas- 
sively advected by the previously prescribed isotropic ho- 
mogeneous and stationary random flow v(r, i), 

^M. = DV^(r, t) - V.[v(r, t)^(r, *)], (44) 

In the standard notation used in ([S|), D stands for the 
"bare" molecular diffusion coefficient. By averaging over 
realizations of v(r, t) we simply obtain from the above 
equation the temporal evolution of < ip{r, t) >. Actually 
taking a (5-like initial condition, this quantity is nothing 
but the probability density for the spatiotemporal dis- 
persion of a unit amount of the randomly advected scalar 
JTqIeoI. Thus the first nonzero moment, 



< rV ><0( r ,t)> = 



dvr l r j < ip(r,t) >, (45) 



R" 



is all that we need to compute an effective diffusion co- 
efficient. The procedure outlined above, although simply 
enunciated, is quite involved in its analytical resolution. 
The best way to proceed is thus to follow the standard 
analytical strategies furnished by the theory of Gaussian 
stochastic processes. The central difficulty arises from 
the non-Markovian nature of the process at hand. Con- 
trolled perturbative schemes are thus necessary. In par- 
ticular a consistent expansion, J2(| based on the smallness 
of the correlation time, to, leads to a closed equation for 
< tp(r,t) >, linear in the autocorrelation tensor, from 
which a diffusive regime is identified through the com- 
mon linear law for the scalar dispersion 

< r 2 > - < r > 2 =< (Ar) 2 >= 4D cS t = A(D + AD)t. 

(46) 

Moreover from such an expansion, the explicit expres- 
sions for the leading contribution to AD can be simply 
evaluated as 



/>oo />oo 

AD C = dsR(0,s) + 4:D dssR"(0, 
Jo Jo 



s), (47) 



where R"(0,s) = d 2 R(r,s)/dr 



Computing these 



left integrals for both spectra we end up with a common 
expression which reads 



AD r 



u t I 1 



v) 



In particular for the Kraichnan's flow 



AD, 



K 



Unt 



o'o 



1 



Dirt 
z 'o 



(48) 



(49) 



An analogous expression would be obtained for the 
Karman-Obukhov's spectrum. 

Both expressions identify the zeroth-order contribu- 
tion, ufao, which can be also viewed as an exactly correct 
limiting case of the classical Roberts analysis for the 
diffusion of a scalar field advected by a rapidly varying 
random velocity field. 



C. Discrete scheme 

Needless to say that given the discrete nature of our 
simulations, the numerical results for the scalar disper- 
sion will be more favorably compared when referring to 
the discrete version of the analytical results given above. 
In particular (to) transforms into 



AD d = 



1 2D 
dsR{A,s)+-^ 



o 



R(2 1/2 A,s) + R(2A,s) -2R( 



For the spectra here analyzed the explicit expression 
reads 



AD d = 



2vN 2 A 4 



E 

fJ-V 



2 ( 27T/A 

V N j 



Q 2 



exp 



2tt/A / 2nv\ 

"*¥ +exp 



2D 



2 A 2 c^ v v 



-1 / 2tt^\ / 2ttv 

_ +cxp f_ l _) +cxp f_ 



D. Results 



Numerical simulation proceeds first by constructing 
the random velocity field with the desired statistical 
properties and then by seeding the system with a 5-like 
initial condition for the dispersed scalar at the center of 
the lattice. Equation ( |i4| ) is numerically simulated by 
using a first-order Euler algorithm for the time variable 
and symmetric forms for the derivative operators. The 



values of A and At where those of Sec. Ill which satisfy 
numeric stability criteria. 



As anticipated in Sec. Ill, the initial condition ?y(r,0) 
is chosen to correspond to the steady state of ?7(r, t) (un- 
derstood here in a statistical sense) (see Eq. (^9|)). In 
this way we are sure to be in an isotropic and homoge- 
neous turbulent environment from the beginning of the 
simulation. Randomly advected by the turbulent veloc- 
ity field the scalar spreads over the lattice, Fig. 2. At 
each time step we measure the variance < (Ar) 2 > and 
we fit, after transients, its temporal evolution to a linear 
law, to obtain D r. 

Summarized in Fig. 3 we present two series of results 
for AD corresponding to two different choices for the 
correlation time of the random flow to. Simulation re- 
sults are compared with the appropriate, perturbative 
obtained theoretical predictions both in their continuous 
and discrete versions. As evidenced in that figure the 
agreement is entirely satisfactory. The general conclusion 
is that a large enough relative correction to molecular dif- 
fusion, up to 25% in Fig. 3, can be accurately predicted 
as long as the correlation time of the random flow, to, 
is small enough relative to the typical time scale for the 



(50) 
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scalar advection Iq/uq. In addition we note that the dif- 
ferences between the discrete and continuous theoretical 
results are easily interpreted when comparing the corre- 
sponding expressions for AD. Actually the correlation 
R(r, s) evaluated at the origin r = or to first neighbors 
may be significantly different given the small value of A 
here employed (A = 1) relative to the elementary spacing 
A (A = 0.5) chosen in all our simulations. 

The role of the parameter to is better analyzed in re- 
lation with the results of Fig. 4. Plotted against the 
inverse of the value of the viscosity, a direct measure of 
t for l fixed as prescribed here, the effective scalar dis- 
persion increases with v~ x . As expected from the intrin- 
sic limitation of the perturbative approach here devel- 
oped, numerical simulations and theoretical predictions, 
although showing similar trends, differ progressively as 
to increases. Complementary to the results depicted in 
Fig. 4, chosen to consider the role of correlation time of 
the random flow, we propose Fig. 5 to examine the influ- 
ence of the correlation length lo- Our theoretical result 
Eq. (4.19) predicts a very small dependence on lo, which 
is the behavior observed in Fig. 5. Note in this respect 
that in going from l = 1.0 to lo = 3.0, eo had to be 
increased nearly two orders of magnitude to keep Uq and 
to fixed, but AD hardly changes. 

The last worthy remark refers to the comparison of K 
and KO spectra. To this end we proposed in Fig. 6 a 
standard representation of AD vs. ufoo for both flows 
with identical integral values lo and to- What we see 
is that the KO spectrum has more statistical dispersion 
and lower values of the effective diffusion than K spec- 
trum. These facts can be understood by comparing both 
spectra in Fig. 1. The most remarkable difference be- 
tween both spectra concerns the inertial subrange they 
mean to represent: KO spectrum is much broader than 
K wavenumber dispersion. This in turn is directly re- 
flected in the different behavior of both spectra for large 
k. Thus KO spectrum shows a richer variety of turbulent 
structures at short distances. For large k, KO spectrum 
still allows structures of the same order than the lattice 
size A. This is not the case of K spectrum where very 
small structures have a very low weight as can be seen in 
the value of the intensity of both spectra for k = 2. So 
in KO spectrum there is less energy concentrated in the 
interval of those k's where the maximum of the energy 
takes place. This would finally result in a reduction of 
the effectiveness of stirring in the case of KO spectrum, 
since small fc-eddies can be thought to be more effective 
in dispersing the scalar. 

V. SUMMARY AND PERSPECTIVES 

An stochastic method, based in a Langevin equation, 
to generate Gaussian synthetic turbulent flows has been 
presented. Characteristic parameters of the flow such as 



its intensity and integral time and space scales are well 
controlled. A relevant aspect of our method is that mul- 
tiple choices of flow spectra can be generated. Two of 
those spectra are explicitly presented here. As a practi- 
cal application of this methodology we have considered 
the study of the diffusion of a passive scalar under the 
influence of two different flow spectra. The role of the dif- 
ferent parameters, in particular the kinematic viscosity, 
which control the characteristic of the flow is discussed. 

The generation of well-controlled flow spectra can be 
a very useful tool in other problems of practical applica- 
tion. For example, this approach was already used in the 
study of phase separation dynamics under stirring [^3| . 
Reactive fronts, i.e., flames propagating under turbulent 
convection is presently under study following also this 
technique 
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FIG. 1. Dynamical evolution of E(k,t) (dashed lines), 
from an initially flat spectrum to its steady state; ti = 0.5to 
and t2 = to; K: Kraichnan's, K-O: Karman-Obukhov. Solid 
lines correspond to the continuous steady spectra, Eqs. (4.1) 
and (4.8), respectively, and dotted lines correspond to discrete 
Eq. (3.11) (uo = 2.5, t = 0.1, Z = 1.5, 100 realizations). 

FIG. 2. Pattern of an initial black drop of the dispersed 
scalar under the influence of the Kraichnan's flow (D — 0.10, 
e = 651.8, A = 2.26, v = 5.1, and t = 160). 



FIG. 3. AD vs u 2 to (Kraichnan's spectrum). Here and 
in the following plots, solid lines correspond to Eq. (4.19), 
broken lines to the discrete expression, Eq. (4.20), and sym- 
bols stand for simulations. Parameter values: D = 0.30, 10 
realizations. Circles: A = 1.00, v — 5.0. Squares: A = 1.69, 
v = 28.5. 



FIG. 4. AD vs the inverse of the kinematic viscosity. 
(u\ = 0.25, A = 2.0, D = 0.30, 10 realizations). 

FIG. 5. AD vs l . {ul = 0.25, t = 0.10, D = 0.30, A = 0.5 
for the open symbols and short dashed line, A = 1.0 for the 
full symbols and long dashed line). 

FIG. 6. AD for Kraichnan (circles) and Karman-Obukhov 
(squares) spectra (to = 0.10, lo = 1.5, 5 realizations). 
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